DVG_Analysis

Kate’s DVG code

message("Loading packages")

Loading packages

library('plyr')
library('tidyverse')
library('reshape2')
library('ggpubr')
library('glue')
library('ggVennDiagram')

# paramters used when running divrge
grouping_param = 5
match_length_param = 28
readLength = 150

# deletion read count cutoffs
count_cut = 30

message("Setting work directory and input file names")

Setting work directory and input file names

wkdir = "/Users/marissaknoll/Desktop/GitHub/Obesity/NewExtractions/H9N2/DVGs"
setwd(wkdir)

if (!dir.exists(glue("{wkdir}/DVG_figures"))) {
        dir.create(glue("{wkdir}/DVG_figures"))
      }

saveitdir = glue("{wkdir}/DVG_figures")

source(glue('{wkdir}/scripts/obese_PlotPrep.R'))

# loading in metadata and coverage data
metafile = glue("{wkdir}/../H9_Metadata.csv")
meta = read.csv(file=metafile,header=T,sep=",",na.strings = c('')) 
meta = filter(meta, resequenced == "yes") %>% filter(cohort != "ND")

transmission_info = "/Users/marissaknoll/Desktop/GitHub/Obesity/NewExtractions/H9N2/TransmissionPairs.csv"
pairs = read.csv(transmission_info, header = T)

meta = merge(meta, pairs, all.x = TRUE) %>% unique()

coverage_passfile = glue('{wkdir}/scripts/H9N2.coverage.pass.check.200.0.95.csv')
cov_check = read.csv(file=coverage_passfile,header=T,sep=",",na.strings = c(''))

# filter for samples that either pass with a yes OR has good average coverage and percentage cov at 200x is > 80
cov_filt_names = cov_check %>% filter(pass == 'YES' | 
                                      mean_coverage >= 200  | 
                                      percentage > 0.4) %>% 
            select(name, segment) %>% 
            unique()

# check segment count
cov_filt_names = cov_filt_names %>% group_by(name) %>% add_tally(name = 'segment_tally') %>% 
                    ungroup() %>% 
                    filter(segment_tally == 8) %>% 
                    unique() 

pull_names = c(levels(factor(cov_filt_names$name)))  # list to pull names from

dvgfile = glue('{wkdir}/H9N2.DVG.FINAL.OneGap.N5.Mis2.M28.G5.csv') # dvg file
dvg = read.csv(file=dvgfile,header=T,sep=",",na.strings = c(''))

dvg = dvg %>% filter(name %in% pull_names) # filter for samples that pass our coverage checks
dvg$sample = dvg$name  # generate new column so we can separate
dvg = dvg %>% separate(sample, c('new','cohort','ferret_id','dpi','rep'), '_')  # separate into info

Warning: Expected 5 pieces. Missing pieces filled with `NA` in 40841 rows [200561, 200562, 200563, 200564, 200565, 200566, 200567, 200568, 200569, 200570, 200571, 200572, 200573, 200574, 200575, 200576, 200577, 200578, 200579, 200580, ...].

CONTROLS = dvg %>% filter(ferret_id == 'HK1073')  # pulling out controls
CONTROLS$rep = CONTROLS$dpi
CONTROLS$dpi = 'stock'  # adding in stock info

dvg = dvg %>% filter(!name %in% c(levels(factor(CONTROLS$name)))) %>% unique()

dvg = rbind(dvg, CONTROLS) # rbind everything so it is all in one dataframe

# prepping rep information
dvg = dvg %>% select(-SegTotalDVG, -new) %>% filter(DVG_freq >= count_cut) %>% unique()  # filter for those that pass cutoffs
rep1 = dvg %>% filter(rep == 'rep1') %>% unique()
rep2 = dvg %>% filter(rep == 'rep2') %>% unique()

# merge reps into one
dvg_reps = merge(rep1, rep2, by = c('cohort','ferret_id','dpi',
                                   'segment','segment_size','strain',
                                   'DVG_group','NewGap',
                                   'NewStart','NewEnd','GroupBoundaries',
                                   'DeletionSize','EstimatedFragLength'), all = TRUE) %>% unique()

# add in zeros
dvg_reps$DVG_freq.x[is.na(dvg_reps$DVG_freq.x)] = 0
dvg_reps$DVG_freq.y[is.na(dvg_reps$DVG_freq.y)] = 0

ggplot(dvg_reps, aes(x=DVG_freq.x, y=DVG_freq.y)) + 
    geom_point() + 
    PlotTheme1

# number of samples?
levels(factor(dvg_reps$ferret_id)) %>% length()

[1] 42

# reorder by segment size
SEGMENTS = c('H9N2_PB2', 'H9N2_PB1',
            'H9N2_PA','H9N2_HA','H9N2_NP', 
            'H9N2_NA','H9N2_MP','H9N2_NS')

cov_check$segment = factor(cov_check$segment, levels = SEGMENTS)

cov_check %>% 
    filter(name %in% pull_names) %>%
    ggplot(., aes(x= segment, y = mean_coverage)) + 
    geom_boxplot() + 
    PlotTheme1

cov_check %>% 
    filter(name %in% pull_names) %>%
    ggplot(., aes(x= segment, y = median_coverage)) + 
    geom_boxplot() + 
    PlotTheme1

cov_check %>% 
    filter(name %in% pull_names) %>%
    ggplot(., aes(x= segment, y = percentage)) + 
    geom_boxplot() + 
    PlotTheme1

df = cov_filt_names %>% select(-segment, -segment_tally) %>% unique()
df$sample = df$name
df = df %>% separate(sample, c('new','cohort','ferret_id','dpi','rep'), '_')

Warning: Expected 5 pieces. Missing pieces filled with `NA` in 10 rows [35, 36, 50, 72, 85, 154, 166, 171, 188, 202].

CONTROLS = df %>% filter(ferret_id == 'HK1073')
CONTROLS$rep = CONTROLS$dpi
CONTROLS$dpi = 'stock'

df = df %>% filter(!name %in% c(levels(factor(CONTROLS$name)))) %>% unique()
df = rbind(df, CONTROLS)

r1 = df %>% filter(rep == 'rep1') %>% select(-new) %>% unique()
r2 = df %>% filter(rep == 'rep2') %>% select(-new) %>% unique()
reps = merge(r1, r2, by = c('cohort','ferret_id','dpi')) %>% unique()

# these are the samples that only had one rep!
setdiff(levels(factor(r1$ferret_id)),
       levels(factor(r2$ferret_id)))

[1] "1411" "1972" "2244"

setdiff(meta$ferretID, reps$ferret_id)  # samples in meta not in seq data

[1] "1411" "1971" "1972" "2235" "2244" "2867" "2868"

setdiff(reps$ferret_id, meta$ferretID) # samples in seq data not in meta

character(0)

m = merge(reps, meta, by.x = c('ferret_id','cohort'), by.y = c("ferretID","cohort"), all = TRUE) %>% 
  filter(inf_route %in% c('Index','Contact','Control'))
         
write.csv(m, glue('{wkdir}/scripts/UPDATED.H9N2.metadata.csv'), row.names = FALSE)

# type check - only stock index direct
print(levels(factor(m$inf_route)))

[1] "Contact" "Control" "Index"  

m$inf_route = factor(m$inf_route, levels = c('Control','Index','Contact'))

m = m %>% filter(name.x != is.na(name.x)) %>% unique()

p1 = m %>% filter(!(is.na(pair_numbers))) %>% unique() %>% 
    ggplot(., aes(x= dpi, y = pair_numbers, fill = diet)) + 
    geom_tile(color = 'black') + 
    PlotTheme3 +
    DietcolScale_fill + 
    facet_grid(pair_diets~inf_route, scales = 'free', space = 'free')

print(p1)

ggsave(p1,
       filename = glue("{wkdir}/DVG_figures/final.samples.pdf"),
       width = 6,
       height = 6, limitsize=FALSE, useDingbats = FALSE)

ggsave(p1,
       filename = glue("{wkdir}/DVG_figures/final.samples.png"),
       width = 6,
       height = 6, limitsize=FALSE) #, useDingbats = FALSE)

dvg_reps = dvg_reps %>% 
            filter(DVG_freq.x > count_cut & DVG_freq.y > count_cut) %>% 
        unique()   # make sure that both reps pass our cutoff

# add in variables for plotting
dvg_reps$ferret_day = paste0(dvg_reps$ferret_id, '_', dvg_reps$dpi)  

m$ferret_day = paste0(m$ferret_id, '_', m$dpi)

stock_temp = dvg_reps %>% filter(dpi == 'stock') %>%
    group_by(ferret_id, cohort, dpi, segment, name.x, name.y) %>%
    add_tally(name = 'seg_deletion_richness') %>%
    unique() %>%
    ungroup() %>% 
    group_by(ferret_id, dpi, name.x, name.y, cohort) %>% 
    add_tally(name = 'deletion_richness') %>%
    ungroup() %>% 
    unique()

s = stock_temp  # will use later

# filter down stock temp information
stock_temp = stock_temp %>% 
            select(ferret_id, dpi, cohort,ferret_day, segment, deletion_richness, seg_deletion_richness) %>%
            unique()


stock_temp = merge(stock_temp, m, by = c('ferret_id', 'dpi','cohort','ferret_day')) %>% 
    unique()

# filter out stock information, calculate dvg richness by segment and across genome for samples
dr = dvg_reps %>% 
            filter(dpi != 'stock') %>% 
            unique() %>% 
            group_by(ferret_id, dpi, segment, name.x, name.y, cohort) %>%
            add_tally(name = 'seg_deletion_richness') %>%
            ungroup() %>% 
            group_by(ferret_id, dpi, name.x, name.y, cohort) %>% 
            add_tally(name = 'deletion_richness') %>%
            ungroup() %>% 
            unique()

# filter down information so you don't have duplicates
richness = dr %>% 
            select(ferret_id, dpi, cohort,ferret_day, segment, deletion_richness, seg_deletion_richness) %>%
            unique()

# merge with metadata info
#richness = merge(richness, m, by = c('ferret_id', 'dpi','cohort','ferret_day'), all.y = TRUE) %>% # CHANGED
richness = merge(richness, m, by = c('ferret_id', 'dpi','cohort','ferret_day'), all.x = TRUE) %>%
    unique() %>% 
  filter(!is.na(inf_route))

# make sure we filter out stock information (will add using the 's' dataframe generated above)
richness = richness %>% filter(dpi != 'stock')

reps_df = rbind(dr, s) %>% unique() # final reps richness df
reps_df = merge(reps_df, m, by = c('ferret_id','dpi','cohort','ferret_day','name.x','rep.x','name.y','rep.y'))  %>% unique() # add metadata

p4 = reps_df %>% 
    select(segment, NewGap, EstimatedFragLength, diet) %>%
    unique() %>%
    ggplot(., aes(x= EstimatedFragLength)) + 
    geom_histogram(color = 'black') + 
    PlotTheme1 +
    labs(x="estimated DVG frag. length (nt)", y='number of unique DVG species') 
print(p4)

ggsave(p4,
       filename = glue("{wkdir}/DVG_figures/deletion.size.pdf"),
       width = 5,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

ggsave(p4,
       filename = glue("{wkdir}/DVG_figures/deletion.size.png"),
       width =5,
       height = 5, limitsize=FALSE) #, useDingbats = FALSE)

p4_alt = reps_df %>% 
    select(segment, NewGap, EstimatedFragLength, diet, inf_route) %>%
    unique() %>%
    ggplot(., aes(x= EstimatedFragLength)) + 
    geom_histogram(color = 'black', binwidth = 50) + 
    facet_grid(inf_route~diet) +
    PlotTheme1 +
    labs(x="estimated DVG frag. length (nt)", y='number of unique DVG species') 
print(p4_alt)
ggsave(p4_alt,
       filename = glue("{wkdir}/DVG_figures/deletion.size.bydiet.bytype.pdf"),
       width = 10,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

lean_index =  reps_df %>% filter(inf_route == 'Index' & diet == 'Lean') %>% 
                unique() %>%
                group_by(NewGap, segment, NewStart, NewEnd) %>%
            add_tally(name = 'lean_deletion_count') %>%
    ungroup() %>%
    select(NewGap, segment, lean_deletion_count) %>%
    unique()


obese_index =  reps_df %>% filter(inf_route == 'Index' & diet == 'Obese') %>% 
                unique() %>%
                group_by(NewGap, segment, NewStart, NewEnd) %>%
            add_tally(name = 'obese_deletion_count') %>%
    ungroup() %>%
    select(NewGap, segment, obese_deletion_count) %>%
    unique()


df = merge(lean_index, obese_index, by = c('NewGap','segment'), all = TRUE) 

head(df)

df$lean_deletion_count[is.na(df$lean_deletion_count)] = 0

df$obese_deletion_count[is.na(df$obese_deletion_count)] = 0

p8 = reps_df %>% filter(inf_route == 'Index') %>% 
                unique() %>%
                group_by(NewGap, segment, NewStart, NewEnd) %>%
            add_tally(name = 'sample_count') %>%
    ungroup() %>%
    select(NewGap, segment,sample_count) %>%
    unique()  %>%
    ggplot(., aes(x=sample_count, y = ..count../sum(..count..))) +
    geom_histogram(color ='black') + 
    labs(x='number of samples with DVG type', y='proportion of DVGs in dataset (index only)') + 
    PlotTheme1 

print(p8)
ggsave(p8,
       filename = glue("{wkdir}/DVG_figures/sample.count.histo.pdf"),
       width = 5,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

ggsave(p8,
       filename = glue("{wkdir}/DVG_figures/sample.count.histo.png"),
       width =5,
       height = 5, limitsize=FALSE) #, useDingbats = FALSE)

reps_df$ave_dvg_freq = (reps_df$DVG_freq.x + reps_df$DVG_freq.y)/2

reps_df = reps_df %>%
  dplyr::arrange(ferret_day, ave_dvg_freq) %>% 
  dplyr::group_by(ferret_day) %>% 
  dplyr::mutate(order_number = row_number()) %>% 
  ungroup() %>%
  unique()

reps_df %>% 
    dplyr::group_by(NewGap, segment, inf_route, diet) %>%
    dplyr::mutate(mean_order = mean(order_number),
          sample_count = n(),
          min_order = min(order_number),
          max_order = max(order_number),
          median_ord = median(order_number)) %>%
    ungroup() %>%
    unique() %>%
    select(segment, NewGap, mean_order, sample_count, min_order, max_order, median_ord, inf_route, diet) %>%
    filter(sample_count > 1) %>%
    unique() %>%
    ggplot(., aes(y=mean_order, x = sample_count)) + 
    geom_point() + 
    PlotTheme1 + 
    facet_grid(.~diet + inf_route)

top_ten = reps_df %>% filter(order_number %in% c(1, 2,3 ,4, 5, 6, 7, 8, 9, 10)) %>% unique()

head(top_ten)

length(levels(factor(top_ten$NewGap)))

[1] 178

max(df$lean_deletion_count)

[1] 112

max(df$obese_deletion_count)

[1] 108

p9 = ggplot(df, aes(x=lean_deletion_count, y = obese_deletion_count)) + 
    geom_jitter(width = 0.1, height = 0.1, alpha = 0.3) + 
    geom_hline(yintercept = 0, linetype = 2, color = 'black') +
    geom_hline(yintercept = 25, linetype = 2, color = 'red') +
    geom_vline(xintercept = 0, linetype = 2, color = 'black') +
    geom_vline(xintercept = 17, linetype = 2, color = 'red') +
    labs(x= 'number of lean samples with DVG', y='number of obese samples with DVG') + 
    PlotTheme1 

print(p9)
ggsave(p9,
       filename = glue("{wkdir}/DVG_figures/sample.count.lean.v.obese.pdf"),
       width = 5,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

ggsave(p9,
       filename = glue("{wkdir}/DVG_figures/sample.count.lean.v.obese.png"),
       width =5,
       height = 5, limitsize=FALSE) #, useDingbats = FALSE)

richness = rbind(richness, stock_temp)
#richness$deletion_richness[is.na(richness$deletion_richness)] = 0

DAYS = c('stock','d02','d04','d06','d08','d10','d12')

richness$dpi = factor(richness$dpi, levels = DAYS)
richness %>% filter(dpi %in% c('d02','d04')) %>%
        filter(inf_route == 'Index' | inf_route == 'Control') %>%
        filter(!(is.na(inf_route))) %>%
    select(ferret_id, dpi, deletion_richness, inf_route, diet, pair_diets, cohort) %>%
    unique() %>%
    group_by(ferret_id) %>%
    add_tally(name = 'n') %>%
    ungroup() %>%
    filter(n >= 2) %>% 
    ungroup() %>%
    unique() %>%
    ggplot(., aes(x=dpi, y = deletion_richness, color = cohort, group=ferret_id, shape = diet)) + 
    #geom_boxplot() + 
    geom_line() + 
    geom_point(size = 2) + 
    PlotTheme1 +
    scale_color_brewer(palette = 'Set1')

    #DietcolScale + 
    #facet_grid(.~cohort)


p7 = richness %>% filter(dpi %in% c('d02','d04','d06')) %>%
        filter(inf_route == 'Index' | inf_route == 'Control') %>%    
    select(ferret_id, dpi, deletion_richness, inf_route, diet, pair_diets, cohort) %>%
    unique() %>%
    group_by(ferret_id) %>%
    add_tally(name = 'n') %>%
    ungroup() %>%
    filter(n >= 2) %>% 
    ungroup() %>%
    unique() %>%
    ggplot(., aes(x=dpi, y = deletion_richness, color = diet, group=ferret_id, shape = diet)) + 
    #geom_boxplot() + 
    geom_line() + 
    geom_point(size = 2) + 
    PlotTheme1 +
    DietcolScale
    #facet_grid(.~cohort)

print(p7)
ggsave(p7,
       filename = glue("{wkdir}/DVG_figures/richness.index.pdf"),
       width = 5,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

ggsave(p7,
       filename = glue("{wkdir}/DVG_figures/richness.index.png"),
       width =5,
       height = 5, limitsize=FALSE) #, useDingbats = FALSE)

colnames(richness)

 [1] "ferret_id"             "dpi"                   "cohort"                "ferret_day"            "segment"              
 [6] "deletion_richness"     "seg_deletion_richness" "name.x"                "rep.x"                 "name.y"               
[11] "rep.y"                 "DPI"                   "sample"                "inf_route"             "diet"                 
[16] "titer"                 "log10_titer"           "HAI_21dpi"             "Ct_Mgene"              "resequenced"          
[21] "pair_numbers"          "pair_diets"           

order_typeday = c('Control_stock',
                  'Index_d02','Index_d04','Index_d06','Index_d08','Index_d10','Index_d12',
                  'Contact_d02','Contact_d04','Contact_d06','Contact_d08','Contact_d10','Contact_d12')

richness$type_day = paste0(richness$inf_route, '_', richness$dpi)
richness$type_day = factor(richness$type_day, levels = order_typeday)

p2 = richness %>% filter(diet == 'Obese' & pair_diets == 'OB>OB') %>%
    select(ferret_id, dpi, deletion_richness, inf_route, diet, pair_numbers, pair_diets, type_day) %>%
    ungroup() %>%
    unique() %>%
    ggplot(., aes(x=type_day, y = deletion_richness, color = pair_numbers, group=pair_numbers)) + 
    #geom_boxplot() + 
    geom_line(size = 1) + 
    geom_point(size = 2) + 
    labs(x='dpi (by index case)', y='DVG richness') + 
    PlotTheme1 +
    scale_color_brewer(palette = 'Set2') #+
    #DietcolScale + 
    #facet_grid(.~inf_route)

print(p2)
ggsave(p2,
       filename = glue("{wkdir}/DVG_figures/obese.to.obese.diversity.pdf"),
       width = 8,
       height = 6, limitsize=FALSE, useDingbats = FALSE)

ggsave(p2,
       filename = glue("{wkdir}/DVG_figures/obese.to.obese.diversity.png"),
       width =8,
       height = 6, limitsize=FALSE) #, useDingbats = FALSE)

gen_rich = richness %>% 
  select(ferret_id, dpi, cohort,deletion_richness, inf_route, diet, pair_numbers, pair_diets, type_day) %>%
  unique()

head(gen_rich)

gen_rich %>% filter(dpi %in% c('d02','d04','d06','stock')) %>%
    filter(inf_route == 'Index' | inf_route == "Control") %>%
    ggplot(., aes(x=diet, y = deletion_richness, group = diet)) + 
    geom_boxplot(outlier.shape = NA) + 
    geom_jitter(width = 0.2, aes(color = diet)) +
    labs(x='segment',y='deletion richness') + 
    PlotTheme1 +
    DietcolScale +
    facet_grid(.~dpi)

Test for significance

o = filter(gen_rich, inf_route == "Index" & dpi == "d06" & diet == "Obese")
l = filter(gen_rich, inf_route == "Index" & dpi == "d06" & diet == "Lean")
t.test(o$deletion_richness,l$deletion_richness)

    Welch Two Sample t-test

data:  o$deletion_richness and l$deletion_richness
t = 0.66564, df = 5.7765, p-value = 0.5313
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -39.95197  69.42816
sample estimates:
mean of x mean of y 
 68.57143  53.83333 

seg_rich = richness %>% #filter(diet == 'obese' & pair_diets == 'obese_obese') %>%
    select(ferret_id, dpi, seg_deletion_richness, inf_route, diet, 
           pair_numbers, pair_diets, type_day, segment) %>%
    unique()

head(seg_rich)

temp = seg_rich %>% select(ferret_id, dpi, inf_route, diet, pair_numbers, type_day, pair_diets) %>% unique()
temp$H9N2_PB2 = 0
temp$H9N2_PB1 = 0
temp$H9N2_PA = 0
temp$H9N2_HA = 0
temp$H9N2_NP = 0
temp$H9N2_NA = 0
temp$H9N2_MP = 0
temp$H9N2_NS = 0

temp = pivot_longer(temp, cols = all_of(SEGMENTS), names_to = 'segment') %>% select(-value)

seg_rich = merge(seg_rich, temp, by = c('ferret_id', 'dpi', 'inf_route', 'diet', 'pair_numbers', 'pair_diets', 'type_day', 'segment'), all = TRUE) 
seg_rich$seg_deletion_richness[is.na(seg_rich$seg_deletion_richness)] = 0

seg_rich$segment = factor(seg_rich$segment, levels = SEGMENTS)
head(seg_rich)

seg_rich %>% filter(ferret_id == 1787) %>% head()

richness %>% filter(ferret_id == 1787) %>% head()

temp %>% filter(ferret_id == 1787)

NA

p3 = seg_rich %>% filter(segment %in% SEGMENTS) %>%
    ggplot(., aes(x=segment, y = seg_deletion_richness)) + 
    geom_boxplot() + 
    labs(x='segment',y='deletion richness') + 
    PlotTheme1 

print(p3)

ggsave(p3,
       filename = glue("{wkdir}/DVG_figures/segment.richness.pdf"),
       width = 5,
       height = 5, limitsize=FALSE, useDingbats = FALSE)

ggsave(p3,
       filename = glue("{wkdir}/DVG_figures/segment.richness.png"),
       width =5,
       height = 5, limitsize=FALSE) #, useDingbats = FALSE)

seg_rich$seg_weight = paste0(seg_rich$segment, '_', seg_rich$diet)
seg_rich$diet = factor(seg_rich$diet, levels = c('Control','Lean','Obese'))

seg_rich %>% filter(segment %in% SEGMENTS) %>%
    drop_na(inf_route) %>% 
    ggplot(., aes(x=segment, y = seg_deletion_richness, group = seg_weight, color = diet)) + 
    geom_boxplot() + 
    labs(x='segment',y='deletion richness') + 
    PlotTheme1 +
    DietcolScale +
    facet_grid(.~inf_route)

p6 = seg_rich %>% filter(segment %in% SEGMENTS & dpi %in% c('d02','d04','d06')) %>%
    filter(inf_route == 'Index' | inf_route == "Control") %>%
    ggplot(., aes(x=segment, y = seg_deletion_richness, group = seg_weight, color = diet)) + 
    geom_boxplot(outlier.shape = NA) + 
    labs(x='segment',y='deletion richness') + 
    PlotTheme1 +
    DietcolScale +
    facet_grid(.~dpi)

print(p6)

ggsave(p6,
       filename = glue("{wkdir}/DVG_figures/segment.index.richness.pdf"),
       width = 8,
       height = 4, limitsize=FALSE, useDingbats = FALSE)

ggsave(p6,
       filename = glue("{wkdir}/DVG_figures/segment.index.richness.png"),
       width =8,
       height = 4, limitsize=FALSE) #, useDingbats = FALSE)

IGNORE THIS SECTION - what it looked like before ferret diets were corrected

# 1973, 1977, and 1986 are listed as obese when they are actually lean, and 1984 is listed as lean when it is actually obese

#seg_rich %>% 
    #filter(!inf_route %in% c('stock','lean','obese')) %>% 
#    drop_na(inf_route) %>%
#    filter(ferret_id == 1408)

#old_obese = filter(seg_rich, ferret_id %in% c(1973, 1977, 1986)) %>% unique()
#old_obese$diet = "Obese"
#old_lean = filter(seg_rich,ferret_id == 1984) %>% unique()
#old_lean$diet = "Lean"
#old_same = filter(seg_rich, ferret_id != c(1973, 1977, 1986, 1984))

#deleteit = rbind(old_obese,old_lean,old_same)

#deleteit %>% filter(!ferret_id %in% c(1973, 1977,1986, 1984) &
#                      segment %in% SEGMENTS & dpi %in% c('d02','d04','d06','stock')) %>%
#    drop_na(inf_route) %>%
#    filter(inf_route == 'Index' | inf_route == "Control") %>%
#    ggplot(., aes(x=segment, y = seg_deletion_richness, group = seg_weight, color = diet)) + 
#    geom_boxplot(outlier.shape = NA) + 
#    labs(x='segment',y='deletion richness') + 
#    PlotTheme1 +
#    DietcolScale +
#    facet_grid(.~dpi)

#seg_rich %>% filter(ferret_id %in% c(1973, 1977,1986, 1984) & dpi == 'd06') %>% 
#  select(ferret_id, segment, dpi, seg_deletion_richness) 

Test for significance

o = filter(seg_rich, inf_route == "Index" & dpi == "d06" & segment == "H9N2_NS" & diet == "Obese")
l = filter(seg_rich, inf_route == "Index" & dpi == "d06" & segment == "H9N2_NS" & diet == "Lean")
t.test(o$seg_deletion_richness,l$seg_deletion_richness)

    Welch Two Sample t-test

data:  o$seg_deletion_richness and l$seg_deletion_richness
t = NaN, df = NaN, p-value = NA
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 NaN NaN
sample estimates:
mean of x mean of y 
        0         0 

# Significance: none

seg_rich %>% filter(inf_route == 'Index' &  segment %in% SEGMENTS) %>%
    ungroup() %>%
    unique() %>%
    ggplot(., aes(x=type_day, y = seg_deletion_richness, color = segment, group=segment)) + 
    #geom_boxplot() + 
    geom_line(size = 1) + 
    geom_point(size = 2) + 
    labs(x='dpi (by index case)', y='DVG richness') + 
    PlotTheme1 +
    scale_color_brewer(palette = 'Set2') +
    #DietcolScale + 
    facet_grid(.~diet + ferret_id, scales = 'free', space = 'free')

p4 = seg_rich %>% filter(diet == 'Obese' & pair_diets == 'OB>OB' & segment %in% SEGMENTS) %>%
    ungroup() %>%
    unique() %>%
    ggplot(., aes(x=type_day, y = seg_deletion_richness, color = segment, group=segment)) + 
    #geom_boxplot() + 
    geom_line(size = 1) + 
    geom_point(size = 2) + 
    labs(x='dpi (by index case)', y='DVG richness') + 
    PlotTheme1 +
    scale_color_brewer(palette = 'Set2') +
    #DietcolScale + 
    facet_grid(.~pair_numbers, scales = 'free', space = 'free')

print(p4)


ggsave(p4,
       filename = glue("{wkdir}/DVG_figures/segment.obese.to.obese.richness.pdf"),
       width = 12,
       height = 4, limitsize=FALSE, useDingbats = FALSE)

ggsave(p4,
       filename = glue("{wkdir}/DVG_figures/segment.obese.to.obese.richness.png"),
       width =12,
       height = 4, limitsize=FALSE) #, useDingbats = FALSE)

End of Kate’s code

Which DVGs are shared between stock and index?

reps_df$DVG = paste0(reps_df$segment,"_",reps_df$DVG_group)

F17_stock = filter(reps_df ,inf_route == "Control", cohort == "F17")
F17_stock_dvg = unique(F17_stock$DVG)
W17_stock = filter(reps_df ,inf_route == "Control", cohort == "W17")
W17_stock_dvg = unique(W17_stock$DVG)
Sm18_stock = filter(reps_df ,inf_route == "Control", cohort == "Sm18")
Sm18_stock_dvg = unique(Sm18_stock$DVG)
Sp19_stock = filter(reps_df ,inf_route == "Control", cohort == "Sp19")
Sp19_stock_dvg = unique(Sp19_stock$DVG)
Sp20_stock = filter(reps_df ,inf_route == "Control", cohort == "Sp20")
Sp20_stock_dvg = unique(Sp20_stock$DVG)

F17_index = filter(reps_df ,inf_route == "Index", cohort == "F17")
F17_index_dvg = unique(F17_index$DVG)
W17_index = filter(reps_df ,inf_route == "Index", cohort == "W17")
W17_index_dvg = unique(W17_index$DVG)
Sm18_index = filter(reps_df ,inf_route == "Index", cohort == "Sm18")
Sm18_index_dvg = unique(Sm18_index$DVG)
Sp19_index = filter(reps_df ,inf_route == "Index", cohort == "Sp19")
Sp19_index_dvg = unique(Sp19_index$DVG)
Sp20_index = filter(reps_df ,inf_route == "Index", cohort == "Sp20")
Sp20_index_dvg = unique(Sp20_index$DVG)

F17_shared = F17_index %>% filter(DVG %in% F17_stock_dvg) %>% filter((DVG %in% F17_index_dvg)) %>% unique()
F17_denovo = F17_index %>% filter((DVG %in% F17_index_dvg)) %>% filter(!(DVG %in% F17_stock_dvg)) %>% unique()

W17_shared = W17_index %>% filter(DVG %in% W17_stock_dvg) %>% filter((DVG %in% W17_index_dvg)) %>% unique()
W17_denovo = W17_index %>% filter((DVG %in% W17_index_dvg)) %>% filter(!(DVG %in% W17_stock_dvg)) %>% unique()

Sm18_shared = Sm18_index %>% filter(DVG %in% Sm18_stock_dvg) %>% filter((DVG %in% Sm18_index_dvg)) %>% unique()
Sm18_denovo = Sm18_index %>% filter((DVG %in% Sm18_index_dvg)) %>% filter(!(DVG %in% Sm18_stock_dvg)) %>% unique()

Sp19_shared = Sp19_index %>% filter(DVG %in% Sp19_stock_dvg) %>% filter((DVG %in% Sp19_index_dvg)) %>% unique()
Sp19_denovo = Sp19_index %>% filter((DVG %in% Sp19_index_dvg)) %>% filter(!(DVG %in% Sp19_stock_dvg)) %>% unique()

Sp20_shared = Sp20_index %>% filter(DVG %in% Sp20_stock_dvg) %>% filter((DVG %in% Sp20_index_dvg)) %>% unique() # still not working
Sp20_denovo = Sp20_index %>% filter((DVG %in% Sp20_index_dvg)) %>% filter(!(DVG %in% Sp20_stock_dvg)) %>% unique() # still not working

stock_shared = rbind(F17_shared,W17_shared,Sm18_shared,Sp19_shared,Sp20_shared)
index_unique = rbind(F17_denovo,W17_denovo,Sm18_denovo,Sp19_denovo,Sp20_denovo)

stock_obese = filter(stock_shared, diet == "Obese") 
o_dvg = unique(stock_obese$DVG)

stock_lean = filter(stock_shared, diet == "Lean") 
l_dvg = unique(stock_lean$DVG)

stock_dvg <- list(Obese = o_dvg, Lean = l_dvg)

StockDVGs = ggVennDiagram(stock_dvg)
print(StockDVGs)
ggsave(StockDVGs, file = "StockDVGs.pdf", path = saveitdir)

Saving 7.29 x 4.51 in image

#ShockSharedDVGs = ggplot(stock_shared, aes(x = dpi, y = DVG)) +
#  geom_point() +
#  geom_line(aes(group = DVG)) +
#  facet_grid(~segment) +
#  PlotTheme1
#print(ShockSharedDVGs)
#ggsave(ShockSharedDVGs, file = "ShockSharedDVGs.pdf", path = saveitdir)

Are there diet-specific DVGs in index ferrets?

index_obese = filter(index_unique, diet == "Obese") 
o_dvg = unique(index_obese$DVG)

index_lean = filter(index_unique, diet == "Lean") 
l_dvg = unique(index_lean$DVG)

diet_dvg <- list(Obese = o_dvg, Lean = l_dvg)

DietUniqueDVGs = ggVennDiagram(diet_dvg)
print(DietUniqueDVGs)
ggsave(DietUniqueDVGs, file = "DietUniqueDVGs.pdf", path = saveitdir)

Saving 7.29 x 4.51 in image

Pulling out diet-specific DVGs

lean = index_lean %>% 
  filter(DVG %in% l_dvg) %>% 
  filter(!(DVG %in% o_dvg)) %>%
  unique()

lean = lean %>%
  group_by(DVG) %>% 
  mutate(count = 1, totalsamp = sum(count))

mult_lean = filter(lean, totalsamp > 1) %>% 
  unique()

obese = index_unique %>% 
  filter((DVG %in% o_dvg)) %>%
  filter(!(DVG %in% l_dvg)) %>% 
  unique()

obese = obese %>%
  group_by(DVG) %>% 
  mutate(count = 1, totalsamp = sum(count))

mult_obese = filter(obese, totalsamp > 1) %>% 
  unique()

lean_uniques = lean %>%
  ungroup() %>% 
  select(segment,DVG_group,GroupBoundaries,totalsamp) %>% 
  unique() %>% 
  arrange(desc(totalsamp))

obese_uniques = obese %>%
  ungroup() %>% 
  select(segment,DVG_group,GroupBoundaries,totalsamp) %>% 
  unique() %>% 
  arrange(desc(totalsamp))

lean_DVG_sizes = lean %>% ungroup %>% select(DVG,DVG_group, DeletionSize, NewStart, NewEnd, diet) %>% 
  unique() %>% 
  arrange(desc(DeletionSize)) %>%
  mutate(name = factor(DVG, levels = unique(DVG)))

#ggplot(lean_DVG_sizes, aes(x = DeletionSize, y = name)) +
#  geom_col() +
  #facet_grid(segment~strain) +
#  PlotTheme1
#ggsave("test_dvg.pdf", a, path = saveitdir, height = 30, width = 5)

#lean_DVG_sizes_plot = ggplot(lean_DVG_sizes, aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
#  geom_segment() +
  #facet_grid(segment~strain) +
#  PlotTheme1  +
#  DietcolScale
#print(lean_DVG_sizes_plot)
#ggsave("lean_DVG_sizes_plot.pdf",lean_DVG_sizes_plot, path = saveitdir, height = 10, width = 10)

lean_DVG_size_seg = lean %>% ungroup %>% select(DVG,DVG_group, DeletionSize, NewStart, NewEnd, segment,diet) %>% 
  unique() %>% 
  arrange(desc(segment),desc(DeletionSize)) %>%
  mutate(name = factor(DVG, levels = unique(DVG)))

lean_DVG_size_seg_plot = ggplot(filter(lean_DVG_size_seg, DeletionSize > 400), 
                                aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
  geom_segment() +
  facet_grid(~segment) +
  PlotTheme1 +
  DietcolScale
print(lean_DVG_size_seg_plot)
ggsave("lean_DVG_size_seg_plot.pdf",lean_DVG_size_seg_plot,path = saveitdir, height = 10, width = 10)

lean_indels_size_seg_plot = ggplot(filter(lean_DVG_size_seg, DeletionSize < 400), 
                                aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
  geom_segment() +
  facet_grid(~segment) +
  PlotTheme1 +
  DietcolScale
print(lean_indels_size_seg_plot)
ggsave("lean_indels_size_seg_plot.pdf",lean_indels_size_seg_plot,path = saveitdir, height = 10, width = 10)

obese_DVG_sizes = obese %>% ungroup %>% select(DVG,DVG_group, DeletionSize, NewStart, NewEnd, diet) %>% 
  unique() %>% 
  arrange(desc(DeletionSize)) %>%
  mutate(name = factor(DVG, levels = unique(DVG)))

#ggplot(obese_DVG_sizes, aes(x = DeletionSize, y = name)) +
#  geom_col() +
  #facet_grid(segment~strain) +
#  PlotTheme1 
#ggsave("test_dvg.pdf", a, path = saveitdir, height = 30, width = 5)

#obese_DVG_sizes_plot = ggplot(obese_DVG_sizes, aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
#  geom_segment() +
  #facet_grid(segment~strain) +
#  PlotTheme1 +
#  DietcolScale
#print(obese_DVG_sizes_plot)
#ggsave("obese_DVG_sizes_plot.pdf",obese_DVG_sizes_plot, path = saveitdir, height = 10, width = 10)

obese_DVG_size_seg = obese %>% ungroup %>% select(DVG,DVG_group, DeletionSize, NewStart, NewEnd, segment,diet) %>% 
  unique() %>% 
  arrange(desc(segment),desc(DeletionSize)) %>%
  mutate(name = factor(DVG, levels = unique(DVG)))

obese_DVG_size_seg_plot = ggplot(filter(obese_DVG_size_seg, DeletionSize > 400), 
                                 aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
  geom_segment() +
  facet_grid(~segment) +
  PlotTheme1  +
  DietcolScale
print(obese_DVG_size_seg_plot)
ggsave("obese_DVG_size_seg_plot.pdf",obese_DVG_size_seg_plot, path = saveitdir, height = 10, width = 10)

obese_indels_size_seg_plot = ggplot(filter(obese_DVG_size_seg, DeletionSize < 400), 
                                 aes(x = NewStart, xend = NewEnd, y = name, yend = name, color = diet)) +
  geom_segment() +
  facet_grid(~segment) +
  PlotTheme1  +
  DietcolScale
print(obese_indels_size_seg_plot)
ggsave("obese_indels_size_seg_plot.pdf",obese_indels_size_seg_plot, path = saveitdir, height = 10, width = 10)

Are DVGs transmitted?

dvg_df = select(reps_df, ferret_day, ferret_id, dpi, diet, inf_route, segment, DVG_group, GroupBoundaries, pair_diets, pair_numbers) %>%
  ungroup() %>%
  unique()
dvg_df$seg_dvg = paste0(dvg_df$segment,"_",dvg_df$DVG_group)

index = filter(dvg_df, inf_route == "Index") %>% unique()
first_time = c("1794_d04","1797_d02","1913_d06","1914_d06","1980_d02","1981_d10","1986_d10","2231_d06","2232_d02","2239_d02")
early_time = c("1794_d04","1797_d02","1980_d02","2232_d02","2239_d02")
direct = filter(dvg_df, ferret_day %in% first_time) %>% unique()
 
samples = unique(index$ferret_day)
  
dvg_transmitted = data.frame()
  
for(i in samples){
  print(i)
  
  n = filter(index, ferret_day == i)
  partner = unique(n$pair_numbers)
  d = filter(direct, pair_numbers %in% partner)
  
  if(nrow(d) > 0){
    
    s = unique(d$ferret_day)
    print(s)
  
    comp = merge(n, d, by = c("pair_diets","pair_numbers","segment","seg_dvg","DVG_group","GroupBoundaries"), all.x = TRUE) %>%
      dplyr::count(pair_numbers,ferret_day.x, ferret_day.y)
    colnames(comp) = c("pair_numbers","Index","Contact","count")
     
    dvg_transmitted = rbind(dvg_transmitted, comp)
   
  }else(print("No transmission"))
  
}

[1] "1408_d02"
[1] "No transmission"
[1] "1408_d04"
[1] "No transmission"
[1] "1408_d06"
[1] "No transmission"
[1] "1412_d08"
[1] "No transmission"
[1] "1414_d04"
[1] "No transmission"
[1] "1414_d06"
[1] "No transmission"
[1] "1416_d02"
[1] "No transmission"
[1] "1416_d04"
[1] "No transmission"
[1] "1789_d02"
[1] "1794_d04"
[1] "1789_d04"
[1] "1794_d04"
[1] "1800_d02"
[1] "No transmission"
[1] "1800_d04"
[1] "No transmission"
[1] "1800_d06"
[1] "No transmission"
[1] "1801_d02"
[1] "1797_d02"
[1] "1801_d04"
[1] "1797_d02"
[1] "1910_d02"
[1] "No transmission"
[1] "1912_d02"
[1] "1914_d06"
[1] "1912_d04"
[1] "1914_d06"
[1] "1912_d06"
[1] "1914_d06"
[1] "1968_d02"
[1] "No transmission"
[1] "1970_d02"
[1] "1981_d10"
[1] "1973_d02"
[1] "1986_d10"
[1] "1973_d06"
[1] "1986_d10"
[1] "1974_d02"
[1] "No transmission"
[1] "1974_d04"
[1] "No transmission"
[1] "1975_d04"
[1] "No transmission"
[1] "1975_d06"
[1] "No transmission"
[1] "1977_d02"
[1] "No transmission"
[1] "1977_d06"
[1] "No transmission"
[1] "1984_d02"
[1] "No transmission"
[1] "1984_d04"
[1] "No transmission"
[1] "2232_d02"
[1] "2232_d02"
[1] "2232_d06"
[1] "2232_d02"
[1] "2232_d08"
[1] "2232_d02"
[1] "2233_d04"
[1] "2232_d02"
[1] "2234_d04"
[1] "No transmission"
[1] "2239_d02"
[1] "2239_d02"
[1] "2239_d06"
[1] "2239_d02"
[1] "2239_d08"
[1] "2239_d02"
[1] "2239_d10"
[1] "2239_d02"
[1] "2240_d04"
[1] "2239_d02"
[1] "2251_d04"
[1] "No transmission"
[1] "2251_d06"
[1] "No transmission"
[1] "2253_d04"
[1] "2231_d06"
[1] "2253_d08"
[1] "2231_d06"
[1] "2254_d02"
[1] "No transmission"
[1] "2254_d04"
[1] "No transmission"
[1] "2254_d06"
[1] "No transmission"
[1] "2254_d08"
[1] "No transmission"

#dvg_transmitted = dvg_transmitted %>% 
#  pivot_wider(names_from = contact, values_from = count)

 contact_dvgs = filter(dvg_df, inf_route == "Index") %>% dplyr::count(ferret_id,dpi,diet) %>%
  ggplot(., aes(x = dpi, y = n, color = diet)) +
  geom_point() +
  geom_line(aes(group = ferret_id)) +
  facet_grid(~ferret_id) +
  ylab("DVG richness") +
  xlab("DPI") +
  PlotTheme1 +
  DietcolScale
print(contact_dvgs)
ggsave(contact_dvgs, file = "contact_dvgs.pdf", path = saveitdir)

Saving 7.29 x 4.51 in image